%% =========================== [Description] ============================
% ---------------- implements the function -----------------
% Under multiple inputs
% Draw the mean variance curve of acceleration displacement         
% Draw a probability density curve for a given second             
% Draw response comparison curves for deterministic vs. random inputs
% ---------------- prefiles -----------------
% quasi0StiffSuspensionCertain.slx
% quasi0StiffSuspension.slx
% n1nDpolyRoots.m
% GetStatistics.m/GetStat.m
% ------------------ enter -------------------
% Cylinder air pressure
% Piston distance from rear wall
% Cylinder node distance
% Effective area of cylinder
%% ###################### [Select excitation signal and output amount] ######################
clear
close all
disp('Pavement sine signal - input 1')
disp('Pavement step signal - input 2')
exc=input('Signal:');
disp('Output select acceleration - input 1')
disp('Output select dynamic deflection - input 2')
disp('Output select dynamic load - input 3')
outputAmount=input('outputAmount:');
swiMat=[-1,1];
excSwitch=swiMat(exc);
%% ############################ [parameter area] #############################
d0Certain=32.5e-3; %% piston distance from back wall m
AcCertain=pi*0.1^2; % cylinder piston area m^2
L0Certain=0.467; % cylinder node distance m
kCertain=3; %cylinder air pressure atm
%--------------------------------------------------------------------------
mCertain=800;
m=800; %Mass on spring kg
m1=100; %Tire mass
Pa=101*10^3;
A_level=16; %-|
B_level=64; %-|-road
C_level=256; %-|
Level=B_level;
v=60; %Vehicle speed
C=2000; %damping
qA=0.1; % excitation frequency
f=1;
%--------------------------------------------------------------------------
svm=4; %random variable sequence number
usvm=[1,2,3]; %deterministic quantity number
D=length(svm); %number of random variables
coeOfVar=[0,0,0,0.1];
MCnumber=500;
vMat=[d0Certain,AcCertain,L0Certain,kCertain];
randomNum=MCnumber^D;
order=3; 
accMCData=zeros(10001,randomNum);
devMCData=zeros(10001,randomNum);
loadMCData=zeros(10001,randomNum);
%% ########################## [compute area] ################################

%% ============================[C] ==================================

sim('quasi0StiffSuspensionCertain.slx')
Ctable={accCertain,devCertain,LoadCertain};
opDataCertain=Ctable{outputAmount};

%% ===========================[PCE] =================================

Collocation=n1nDpolyRoots(order,D); 
groupNumber=length(Collocation(:,1));
for groNum=0:groupNumber  
    if groNum==0
        [d0,Ac,L0,k]=iniRanVar(d0Certain,AcCertain,L0Certain,kCertain,...
                     coeOfVar,Collocation,groupNumber,1);
        sim('quasi0StiffSuspension.slx')
        outputData=zeros(length(a.time),groupNumber);
    else
        [d0,Ac,L0,k]=iniRanVar(d0Certain,AcCertain,L0Certain,kCertain,...
                     coeOfVar,Collocation,groupNumber,groNum);
        sim('quasi0StiffSuspension.slx')
        switch outputAmount %Select the output amount
            case 1 
                output=a.data;
                 yLabelText='\it{Acceleration}$\rm{(m/s^2)}$';
                proUnit='$\rm{(s^2/m)}$';
            case 2
                output=fd.data;
                yLabelText='\it{Deflection}$\rm{(m)}$';
               proUnit='$\rm{(1/m)}$';
            case 3
                output=Fd.data;
                yLabelText='\it{Relative Dynamic Load}';
                proUnit='';
            otherwise
                error('Select 1\2\3');
        end
        outputData(:,groNum)=output;
    end
   %-----------------------------------------------------------------------
    proPers=floor(groNum/groupNumber*1); %
    clc %---- progress bar area
    disp(['PCE[','#'*ones(1,proPers),' '*ones(1,100-proPers) ...          %
        ,']',num2str(proPers),'%']) %
   % -----------------------------------------------------------------------
end
pointNumber=200;
timeNumber=10;
sampleTime=length(a.time)-1;
stepl=timeNumber/pointNumber;
time=0:stepl:timeNumber-stepl;
PCmean=zeros(pointNumber,1);
PCvar=zeros(pointNumber,1);
MCdata=zeros(pointNumber,MCnumber^D);
MCinput=randn(MCnumber,D);
for i=1:pointNumber    
    [PCmean(i),PCvar(i),MCdata(i,:)]=PCE_GS ...
    (Collocation,outputData((sampleTime/pointNumber*(i-1))+1,:)',order,MCinput);
% [PCmean(i),PCvar(i)]=GetStat ...
% (Collocation,outputData((sampleTime/pointNumber*(i-1))+1,:)',order); 
    %----------------------------------------------------------------------
    proPers=floor(i/pointNumber*99); %
    clc %---- progress bar area
    disp(['PCE[','#'*ones(1,1),'#'*ones(1,proPers),' '*...                %
        ones(1,99-proPers) ,']',num2str(1+proPers),'%']) %
    % ----------------------------------------------------------------------
end

%% ============================ [MC] =================================

for index=1:randomNum
    [d0,Ac,L0,k]=MC_initialize(vMat,coeOfVar,svm,...
    index,MCnumber,MCinput);
    sim('quasi0StiffSuspension.slx')
    accMCData(:,index)=a.data;
    devMCData(:,index)=fd.data;
    loadMCData(:,index)=Fd.data;
    %----------------------------------------------------------------------
    proPers=floor(index/randomNum*100); %
    clc %---progress bar
    disp(['MC [','*'*ones(1,proPers-1),'>',' '*ones(1,100-proPers) ...    %
        ,']',num2str(proPers),'%']) %
    % ----------------------------------------------------------------------
end
MCtime=a.time;
MCtable={accMCData,devMCData,loadMCData};
opDataMC=MCtable{outputAmount};
MCmean=mean(opDataMC,2);
MCsigma=sqrt(var((opDataMC')));
MCsigma=MCsigma';

%% ########################## [drawing area] ############################### 

%% ====================== [standard deviation and mean plot] ============================
 close all
figure(1)
errorbar(time(1:2:199),PCmean(1:2:199),sqrt(PCvar(1:2:199)),'-s',...
    'color',[0,0,0],...
    'linewidth',0.7,...                        % Set line width     
    'MarkerSize',5,...                         % Set marker symbol size
    'MarkerEdgeColor','black',...                % Set marker edge size
    'MarkerFaceColor','red') %set the marker internal size
%axis([0,10,-30,30])
xlabel('Time(s)')
ylabel(yLabelText)
disp('Program operation finished')
%% ==================== [Chaotic mean and model mean plot lines] ======================
close all;
axisMax=1.5*max(PCmean);
figure(2)
plot(opDataCertain.time,opDataCertain.data,'LineWidth',1.5,'color','red')
hold on
plot(MCtime,MCmean,'-.' ,'MarkerSize',1,'LineWidth',1.5,'color','blue')
hold on
plot(time,PCmean,'x','MarkerSize',9,'LineWidth',0.7,'color','black')
ylim([-0.7*axisMax,axisMax])
ylabel(yLabelText,'Interpreter','latex','FontSize',20)
xlabel('Time(s)','Interpreter','latex','FontSize',20)
set(gca, 'Fontname', 'Times New Roman','FontSize',20);
legend('Deterministic Input','Random Input(MC)','Random Input(PCE)')
%% ======================== [probability distribution for a given second] ============================
tm=[17 59 83]*0.1;
for l=1:length(tm)
    t=tm(l);
    tNumber=floor(t/stepl+1);
    dth=50*tNumber-49;
    [PCEpdf,PCEt] = ksdensity(MCdata(tNumber,:));
    [MCpdf,MCt] = ksdensity(opDataMC(dth,:));
    figure(3)
    subplot(1,3,l)
    plot(MCt,MCpdf,'LineWidth',3,'color',[30,144,255]/256)
    fill(MCt,MCpdf,[192,192,192]/256)
    hold on
    plot(PCEt,PCEpdf,'LineWidth',1,'color',[255,69,0]/256)
    xlabel([yLabelText,'-', num2str(tm(l)),'s'],'Interpreter','latex','FontSize',20)
    ylabel(['Probability Density',proUnit],'Interpreter','latex','FontSize',20)
    legend('MC','PCE-MC')
    set(gca, 'Fontname', 'Times New Roman','FontSize',20);
end
%% ================= [Mean and random input probability density] =========================
tm=[22 27 18]*0.1;
for l=1:length(tm)
    t=tm(l);
    tNumber=floor(t/stepl+1);
    dth=50*tNumber-49;
    [PCEpdf,PCEt] = ksdensity(MCdata(tNumber,:));
    [MCpdf,MCt] = ksdensity(opDataMC(dth,:));
    figure(4)
    subplot(1,3,l)
    plot(MCt,MCpdf,'LineWidth',3,'color',[30,144,255]/256)
    fill(MCt,MCpdf,[192,192,192]/256)
    hold on
    plot(PCEt,PCEpdf,'LineWidth',1,'color',[255,69,0]/256)
    hold on
    pdfx=opDataCertain.data(floor(t*1000))*ones(2,1);
    pdfy=[min(PCEpdf),max(PCEpdf)];
    plot(pdfx,pdfy)
  
    xlabel([yLabelText,'-', num2str(tm(l)),'s'],'Interpreter','latex','FontSize', 20)
    ylabel(['Probability Density',proUnit],'Interpreter','latex','FontSize',20)
    legend('MC','PCE-MC')
    set(gca, 'Fontname', 'Times New Roman','FontSize',20);
end    















